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ABSTRACT 


A numerical model was developed depicting the wind-driven 
circulation in an ocean basin containing an island. This 
linear, Е filtered model was utilized to test and 
evaluate the "Hole Relaxation" technique (Allen 1954) in 
preparation for future comparative studies with a free sur- 
face model and later incorporation in a multi-level world- 
ocean model with islands. Eight different data cases were 
studied to evaluate the model's ability to properly treat a 
variety of island sizes and locations. It was found that as 
the relative size of the island was increased at a specific 
location in the western half of the basin, the volume trans- 
port per unit depth between the mainland and the island, as 
well as the total kinetic energy at the end of a 210 day in- 
tegration decreased; if the island size was decreased, the 


parameters' magnitudes were noted to be larger. 
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I. INTRODUCTION 


For decades, men driven by the reality of the world's 
oceans and their geometric configurations have strived to 
model the physical boundary characteristics and the processes 
which govern the ocean-atmosphere realm. These men have 
recognized the future needs of the world's forecasting media 
by developing models which have answered some old queries 
and have further related new theories to make present-day 
phenomena more easily understood. Most studies were con- 
Siderably enhanced by the advent and sophistication of 
numerical techniques. 

In 1947, Sverdrup showed for the first time the major 
role of the curl of the wind stress in determining the 
general mid-ocean distribution of the meridional current. 
Later, Stommel (1948) demonstrated the decisive role of the 
latitudinal change in the coriolis parameter in the forma- 
tion of intensive boundary currents along the eastern shores 
of continents. Stommel also used the linearized vorticity 
equation with linear bottom friction. In 1950, Munk, utiliz2- 
ing lateral frictional dissipation, directed his studies at 
more detailed features of linear viscous boundary currents, 
all the while restricting these studies to simple domains of 
integration (i.e., rectangle and triangle) and by zonal wind 
fields. All of these models used a vorticity equation con- 


 Sidered on a Beta plane. 


2 





~“Purther studies, now advanced by numerical techniques, 
were undertaken. Gates in 1968 showed the development of 
the large-scale wind-driven circulation using the primitive 
equations of motions for a homogeneous ocean model charac- 
terized by a free upper surface. This time-dependent model 
permitted not only easy observation of the individual tran- 
sients' role in the circulation development, but also exten- 
sive examination of the role of divergence in the large-scale 
circulation by the use of the free surface upper boundary 
condition. Bryan in 1968 and later in 1969 with Manabe, by 
means of a joint atmosphere-ocean numerical model with a 
rigid upper boundary in the ocean, investigated the role of 
the ocean circulation in producing realistic features of the 
climate of the atmosphere. 

Also in 1969, Takano studied the barotropic ocean circu- 
lation in a doubly connected world ocean using the vorticity 
equation (non-stationary) and the numerical method of "hole 
relaxation."  Kamenkovich et al. (1969), in an independent 
effort, further displayed successful simulation of actual 
shoreline profiles by numerical methods in a complete circu- 
lation of the world ocean (stationary). In 1971, Alexander 
extended Gates' work to a baroclinic free surface model depict- 
ing a doubly connected domain, with islands in the western 
boundary region. 

The primary purpose of the following study is to test a 
numerical technique for calculating the large-scale wind- 


driven circulation in an ocean basin which contains an island. 





It Is not the purpose of this model to physically depict any 
particular ocean or seasonal wind pattern; rather, the goal 
15 to simplify the circulation problem in order to study and 
evaluate the numerical technigue at the sacrifice of precise 
world simulation. 

The numerical solution of a barotropic ocean with an 
island has been achieved in both linear and non-linear models. 
The present work is totally a linear effort utilizing the 
vorticity eguation which filters the external high freguency 
gravity waves. This permits the use of a longer time step 
(<14 hours) than in the free surface (non-filtered) models. 
The application of the "Hole Relaxation" technique (Allen 
1954) with respect to the island boundaries (interior 
boundaries to the ocean model) is employed with a major goal 
to test and evaluate the techniques on this simple model for 
future development. Furthermore, the results of this 
evaluation will be used in calculating the barotropic mode 
in more sophisticated world ocean models which can ill- 
afford the luxury of short time steps as required in the 
free surface models. 

The present study did not attempt to argue in favor of 
a filtered model rather than a free surface model for the 
calculation of the barotropic mean currents. Its goal is to 
simply test a technique whereby the merits of the two 
numerical approaches (free surface vs. rigid surface) may be 


compared. 
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ТТ. THE MATHEMATICAL STATEMENT OF THE PROBLEM 


A. VORTICITY EQUATION DERIVATION 

The problem is to obtain a solution to the wind driven 
ocean circulation in a rectangular ocean basin of uniform 
depth containing an island. 

The linearized equations of horizontal motion (1) and 
(2) together with the continuity equation (3) for a homo- 


geneous fluid are as follows: 

















дч _ 1 др _ 2 1. х 
ae fv) + D n AV^u + um (1) 
о о 
dv 1 ie ab 
ара + SR = ду2у + — (2) 
dt p ду 202022 
O о 
ди ду ðw _ 
3x T dy ал am er 0 (3) 


Friction forces are crudely represented by a lateral 
eddy diffusion of momentum (AV^u, AV?v) and a vertical eddy 
stress Ge): All other symbols have their usual meaning. 
Cross differentiation of (1) and (2) to eliminate pressure 


results in the following vorticity equation: 


д ,9У du du О 
3t ax й 3 y? Ш> Шат Б Sy? E 
AT AT 
9 y2, -2 y? 19 _у__х 
А(5- усу ду V ust р. 3, (52 ду , (4) 
where ß = Sn constant 
ду 


By vertically integrating the continuity equation (3) from 
the bottom of the ocean (z = -D) to the mean sea 1еуе1 (2-0), 


we have 
w(-D) 


+ Ој 
7 


> 0 (5) 


а2 | а> 
м [с 1 
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where (u,v) is the yertically integrated velocity defined by 


о 
и = = Ее” 2s 
-D 
and 
= 197 
у = т” v dz . (6) 
-D 


The boundary conditions on w(-D) and w(o) are 


0 


w(-D) 


м (о) 


027. (7) 

The first boundary condition in (7) is kinematic and 
applies at the flat ocean floor. The second condition is 
the approximation which removes external gravity waves from 
the model and thereby allows longer time steps to be 
utilized. The inaccuracy of the solution due to the balance 
approximation of requiring the vertical velocity to be zero 
at the sea surface can be examined by comparison of the 
results of this present calculation with similar calculations 
using a barotropic primitive equation model with a free sur- 
face as used for example by Gates in 1967. As a result of 


(7), (5) becomes 

du , dv o | 

9x y dy Шо c8) 
Consequently, a streamfunction (x,y) is defined for the 


vertically integrated velocity as follows: 


0 

u Jy 

- 5 | 

y= 2 (9) 


The vertical integration of the vorticity equation (4) апа 


the use of (8) and (9) results in 
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3 729 68-33 = дт? (угу) 


9x 
OT OT IT OT 
ЕС - n 0-Х - =) ы 
o D' 5x dy “2-07 Эх ду "z--D. ' 


(10) 


Ihe boundary conditions on the stress at the ocean floor 


(2 = -D) and at the sea surface (z = 0) are as follows: 
At z = -=D 
Тік E = 0 (11) 
while at z = 0 
== ШУ; 
te F cos G ) 
Б = 0 , (12) 


where F is the amplitude of the zonal components of the wind 
Stress and b is the north-south extent of the ocean basin. 
This model in the islandless case was a time dependent repre- 
sentation of Munk's model (1950). Therefore, the results 

can easily be compared with Munk's analytic solution as a 
check on the accuracy of the numerical solution. 

Thus, there is no bottom stress, and the surface stress 
corresponds to a pattern of westerly winds in the northern 
half of the basin, and to easterly winds in the southern 
half of the basin. In addition, it is important to have a 
stress that has a non-zero curl since only the curl of the 
wind stress actually does the forcing in the vorticity equa- 
tion. Using (11) and (12), the vorticity equation takes the 
final form, 

у250 = pl - a sin EL) £AV?(V?y) . (13) 
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В. BASIN COASTLINE AND ISLAND BOUNDARY CONDITIONS 

The boundary conditions along the island and basin coast- 
lines (mainland) extremities used in the solution of (13) are 
those of zero normal flow and of zero slip. 

The condition of no normal flow бул = 0) was implemented 
Бу requiring that the streamfunction be equal to zero on the 
ocean periphery and that the same function be a constant 
function of time on the island's perimeter. These conditions 


ате written as: 


O (mainland boundary) 


ф 
Y 


c(t) (island boundary) (14) 

The value of the constant c(t) determines the net flow 
(volume transport) between the mainland and the island and 

is calculated by the "hole relaxation" procedure (Allen 1954) 
as a part of the solution. 

The conditions of zero slip and no normal flow are 
implemented in the calculation of the lateral friction force 
term on the right hand side of (13). These conditions are 
"written as: 


u, v = 0 (all island and mainland 
boundaries). (15) 


Details of the formulation of these boundary conditions are 
contained in section three where the finite difference equa- 


tions are given. 


С. PHYSICAL CHARACTERISTICS OF THE NUMERICAL MODEL 
The model presented here is a rectangular basin of 


length, L = 8000 km; of breadth, b = 4000 km; and of constant 
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depth, Depth = 2 km. (See Figure 1.) The total grid size is 
4l x 21 in the horizontal 25 the computational grid size 
is 40 x 20. The distance between adjacent grid points is 
200 km in both x and y. The value chosen for Beta makes the 
grid centered at 32.5°N. Figure (1) also displays a portion 
of the 861 nodal points and enous how an arbitrary island 
May appear. The arbitrarily sized, rectangularly shaped 
island may be placed at any selected location throughout the 
grid. The model is not strictly designed to simulate any 
particular ocean or ocean currents system, but rather is 
utilized to evaluate the solution techniques to be tested 
and secondarily to evaluate the evolution of the physical 


flow patterns around the island with respect to time. 


D. EXPERIMENTS CONDUCTED 

Initially, an islandless ocean model was studied. This 
model utilized the same equation (13) as its basis of solu- 
tion. Again, this model was analogous to the solution of 
Munk's model (1950) when the numerical solution achieved 
steady-state conditions. All characteristics were as de- 
scribed above with the exception of an island. A duration 
of 210 days was selected as the run time for this phase in 
order to be able to compare all possible "island" integra- 
tions with this islandless case. 

Next, an island-ocean basin phase was introduced by 
adapting the boundary conditions and restrictions of the 


hole relaxation technique to be tested (Allen 1954). A 


о 
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time plot of the total kinetic energy showed that 210 days 
was adequate to reach a statistical steady state, and to 
give an accurate comparison to the steady-state characteris- 


tics of the islandless phase results. 
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ШЕПТЕ DEP PERENCE EQUATIONS 


Individual terms on the right hand side of (13) are ex- 
pressed in finite difference form in the following paragraphs. 
The grid points in the annulus are discussed first in each 
case, then follows a discussion of the grid points on the 
island perimeter. Island corner points and the points on the 
perimeter adjacent to the corners are not discussed. Readers 
are referred to the documented computer program in Appendix A 
for further details. It should also be clarified that all 
terms on the right hand side of (13) were calculated after 
multiplying the entire equation by the area represented bya 
grid point, namely (Ду The reason for this action was 
that the time rate of change of vorticity multiplied by the 
area represents the time rate of change of the circulation. 
This time rate of change of the circulation was the variable 
"F1" as seen in Appendix A. 

It should be pointed out that a forcing function of Fl=0 
was assumed (Allen 1954) in the interior nodes of the island. 
Finally, since it was assumed that Ах = Ау throughout the 
grid, the term "Ax" was utilized throughout the subsequent 


discussion and all sections of this thesis. 


А. THE BETA TERM 
In the annular locations while maintaining boundary con- 
dition (14), the beta term was expressed as 


Ут+1 7 ?т-1 


-8 30 ax? => -B ( ) Ах (16) 
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{i the beta terms on the north and south perimeters were 
considered to be zero due to boundary condition (14), namely, 
zero normal flow, then the boundary conditions of constant 
values of JU on the perimeter of the island and zero values of 
UV on the mainland, inherent in (14), had to be maintained on 
the east and west island perimeter nodal points. The con- 
Straint that the integrated sum of all beta terms along any 
latitudinal plane or axis was to be zero, guided the formula- 
tion of the finite difference form of the beta term on the 
east-west perimeters. This integral constraint which was a 
property of the differential expression for got, also had to 
hold true in the finite difference analogues, as long as the 
meridional velocity component on the island perimeter was 
evaluated by means of an uncentered difference approximation. 
Therefore, when the beta term was calculated on the island 
perimeter as described above, the uncentered finite difference 


forms used were: 

оош Z. Е 

В Ах epu. m WI Ах , (17) 

for the western perimeter, and 

20007 aa - 
for the eastern perimeter. 
B. CURL OF WIND STRESS TERM 

In both the case of the annular location and that of the 


island perimeter locations, the wind stress term was defined 


in the same manner, namely: 


LE 











- EX б сіп Geo => - re Ae sin (Y) , (T9) 
о O 
у = cil (20) 


C. CURL OF FRICTION FORCE TERM 

Boundary conditions utilizing zero slip and zero normal 
flow (14) were again selected to enable dissipation to occur 
at the boundaries through the linear curl of the friction 
force term in the model. 


The friction force term was expressed as follows: 
2 2 > 2 € ЩЕ > 2 9 y2 
Е 221, 7,701: 021) 


The finite difference form in the open oceanic areas 


away from the boundaries was based on (21) as 


У? у +V?v И TV^v 
2 LJ Ё 1-1 ПЕ I-1,J-1 
[ CAV? c) NL АИ —— ee) 
LR] 2 2 
У? и +V7u V^u +7%u 
MJ I-l,J Edel I-1,J-1 
AA (AAA DP (22) 
where 
2 Ут+1 ШЫ Ут+1 N о ru 1+1 тај 
V^v CIEL S2 7 pl таз лал V, 11 7 "T1 I1, , 
IS 
(22a) 
with a similar expression for V^u » and where u and 
[ST IRE 
MSS are defined in terms of Ше аз, 
Loa Ах 2 2 : 
ап4 
ty 
НЕ ОСЕ И тет о 210-1 257 
ELE iC => =) =( : ТС” 
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On the island perimeter, the same basic finite difference 


form (21) applied except that the following alterations 


existed: 


i) on the northern and southern perimeters, 
DM (23) 
and ii) оп the eastern and western perimeters, 
I sago (24) 
ду 
These conditions preserved boundary condition (14) with 
respect to no normal flow. 
Having made use of (23) and (24) along with the zero slip 
condition, the following finite difference forms resulted 


for use in (22) on the island perimeter: 


i) southern perimeter of the island 


(22)% 5 [-ачу ј- • 


1 о u + u t 


тг. 
(25) 


іі) northern perimeter of the island 


A . 
Lac uM Ur Tere T1  T«2 943 11, 42^ 142, 9627 ! 
(26) 


iii) western perimeter of the island 


l À ! 
— га dm => а= 9 2 
сва) mu E A Desi Ут,у+2 йг! tei, J Ут, 1-11} (27) 


and iv) eastern perimeter of the island 


A 
A 21217/7:11,5:2771:2,141"“142,3"“141,3-17" 


(28) 
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This complex perimeter calculation has been clearly 
illustrated in Figure 2. 

Suppose the point (1,J) in Figure 2 was assumed to be 
located on the western perimeter of the island. The rela- 
tionship, for example, 

2 2 2 2 
Ži, , " TV, - E (29) 
demonstrated how the staggered grid points, numbered one 
through 12, were used to compute б у : 

The right hand side of (29) was written as the sum 
2 


1,3 
factor assigned to v at the nodal point (I + i, J + j) in 


ТН НІН where т was the net weighting 


the computation of 2 үгу at the nodal point (I,J) in the 
island perimeter. Having rewritten the right hand side of 


(29), and having used the zero slip condition, шаг 


Maši gi; it was found that: 


(30) 


V - V 


ШАП СУТ utt MRI ges 


eS ТЕТ ІІІ. 
Readers should note that in the staggered grid depicted, 


was located at point number one, was located at 


Ут, 1+1 


point number two, etc., and that the terms of (30) were 


ТЈ 


signed as a result of (29),i.e., points one and two had 
w = +4, while points nine through 12 had w = -l. 

The above served to demonstrate the derivation Of (27: 
Equations (25), (26), and (28) were derived in a similar 


fashion. 
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Island Perimeter 


Figure 23 


Key: 
© = Grid Nodal Points 
Staggered Nodal Points 


V 


2 2 
O ими, уу 


о 
" 


Illustration represents a portion of 
іле co putational orid, wich includes 
the island perimeter. The side of the 
grid which represents ocean or island 
depends on the case to be evaluated, 
Че о ШО БТ Lor the Equation 
(29) cited. The shifting and/or re- 
orientation of the island portion of 
the figure will permit north, south, 
east or west perimeter calculations. 
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D. OTHER TERMS 

In addition to the three terms on the right hand side of 
(13), the finite difference forms for the residual calcula- 
tion in the relaxation phase of the model's solution were to 
have been expressed. Yet prior to that discussion, a brief 
description of the time differencing procedures was formu- 
lated. 

For accuracy, a centered time difference scheme was used. 
This was the so-called "leapfrog" scheme which was charac- 
terized by the desirable property of being neutral and the 
undesirable property of producing a computational mode in 
time for the system solution. This resulting mode, caused 
by the presence and utilization of three time levels in the 
description of m , was removed by the periodic employment 
of a Euler-Backward (Matsuno) time scheme. The Euler-Back- 
ward scheme was used initially at the first time step and 
again at every 50 time steps thereafter to prevent separation 
of the solution in time. This scheme removed the computa- 
tional mode by utilization of only two time levels for the 
calculation of oe 

The relaxation phase of the model solution had an equa- 


tion that was to be solved which was expressed as, 


Wig (2%) - Е1 = 0 (31) 


where Fl was expressed by the right hand side of (13) multi- 
plied by О ТЕ it is assumed that the symbol Y' was the 
first derivative of Y with respect to time, then the relaxa- 


tion residual, RESID, was expressed as 
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CES ! ! ! -41)! - 
СКЕБІр)г у= Фуру th pea tha, seater, a ti, dell, zo 


(32) 
in both the ocean basin grid points and at the island grid 


points. 


E. RELAXATION DESCRIPTION 
1. Open Ocean Area 
The Liebmann over-relaxation technique was utilized 
in all open ocean areas of the computational grid. This 
technique was found to have the most. fee aera application 
to the model in which forcing by large scale stress dominated. 
Following the standard relaxation technique, if it 


ОООО the solution of (31) the absolute value of 


RESID, 1 іп (32) was greater than some small allowable error, 
3 
E 
€, a new Y! denoted by vy! was defined as follows: 
PES LI 
Peann E cndi (33) 
[ST 1,3 T5 7 
where the correction С. j was determined by the condition 
* 
that the new residual value RESID, I calculated using Ur 1? 
3 > 
must be a negative fraction of RESID, I? im en 
3 
* 2 
= - : 34 

RESIDI y о RESID, , (34) 
A value of a = 0 therefore corresponds to no over- 
relaxation. Since 

* * 
= ' ' ' ! = ! -Е1 35 

DESIDI J Er mp qe E a цэг С ES (35) 


is true, then it was easily found that the relationship for 


the correction C was 
ТЕМ 
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го. и RESID, | (36) 

An optimum value for the over-relaxation coefficient, 
a, of 0.75 was selected for use in these open oceanic regions. 
This value was the result of several test evaluation runs 
and was rather large due to the large scale nature of the 
forcing and resultant motion. More details as to the appli- 
cation of this technique to the open oceanic areas of the 
grid can be found in the ample commentary of the computer 
program located in Appendix A. 

2. Island Areas 

In the island areas, the "Hole Relaxation" technique 
was used (Allen 1954). This technique assumed that the same 
Equation (31) was to be solved with the exception that the 


island points were cumulatively treated like a block where 


the residual on the island, RESIDI, was defined as 


RESIDI у. (vV?^y'- Fly y) : (37) 


island 


! ! 


where the symbol represented a four grid point 


average, and the symbol ' | ' represents the finite dif- 
i island 


ference approximations to an integral over the island of 
the individual residual values. This total residual value, 
RESIDI, was made to approach zero in compliance with the 
hole relaxation technigue in order to achieve an accurate 
solution for (31). As a result of the averaging process 
used, RESIDI was written as the Әне of the individual nodal 


point residuals, each multiplied by a weighting factor where 
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the corner perimeter nodal points were assigned a factor of 
0.25, while the lateral perimeter points were given a factor 
of 0.50, and the interior points a factor of 1,0. 

If the subscripts 1,2,3,4; E,W,N,S; and INT. were 
used to denote the contribution to RESIDI from the four 
corners, the east, west, north and south lateral perimeter 
points and the interior nodal points respectively, then 
since Fl < 0 on the 20 interior, (37) was expanded as: 


RESIDI-0.25[(V^Q'-F1),*(V^V'-F1),*(V^y'F1) ,F (V^ j- F1) , ] 
*£0.50[Z(V?^y'-rF1)«-Z(V?^y'-F1)4Z(V^g'-F1)4Z(V?j'-F1)] 
E W N S 


+1.0[ 2 (у?ф')] . (38) 
INT. 


In accordance with the hole relaxation technique 
(Allen 1954), and the boundary conditions of no normal flow 
(14), ф' had the same value at all points on the island, i.e., 

I,J = p' (on the island) . (39) 

The object of the hole relaxation technique was to 
find the value of wW' for which the absolute value of RESIDI 
was less than the error tolerance e . The procedure used 
to do this was as follows. After completing each relaxation 
sweep over all interior ocean points, RESIDI was calculated 
from (38) and its absolute value tested. If it was equal to 
or less than € , no correction was added to Y', and if no 
correction had been made to any of the interior ocean points 
on that last sweep, then т, and y' comprised the solution. 
If the absolute value of RESIDI was larger than € , a correc- 


tion, CI, was added to y' and a new relaxation sweep over 
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all grid points was initiated. The correction CI was deter- 
mined in such a way that the new RESIDI would be a negative 
fraction of the current RESIDI. The correction CI for the 
island block was determined similarly to the development of 
EL J for the ocean points as follows: 

It was first noted that in (38), 


a. for each corner grid point, 


(v?y') = -2y'+ non-island Q7 , ; (40) 
n=1,2,3,4 2 


b. for each group of lateral nodes not on a corner, 
for example the eastern perimeter, 


E(QV^ y") => (-р')+»; (поп-151апа эг J? (41) 
E E E 1 


L(V5y')2(-9') (N0,) *Z (non-island O d (42) 
E E , 


where NO, was the number of non-corner nodal points on the 
island's eastern perimeter. Finally 
с. for the nodal points on the interior of the 


Es land, 


(Wy, yo. (43) 


If the domain of the island was described by desig- 
nated I,J indices such that the minimum and maximum values 
of I and J were IMIN and IMAX, JMIN and JMAX respectively, 
then the number of nodal points on a particular island side 
were expressed in terms of these indices. The number of 
non-corner lateral nodal ee the example stated in 


(42) was, 
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NO, z JMAX - (JMIN t1 1) . (44) 


Similar expressions for westerly, northerly, and 
southerly lateral side descriptions were evolved for МО у, 


NO and NO respectively. 


В > ©? 
When Equations (40), (42), (43), and (44) and their 
associated counterparts for the other perimeters were used 


to expand (38), the resulting expression for the island 


residual was, 


RESIDI = 0.25[-8y' + x (non-island фт pa? (Fl) ] 
n=1,2,3,4 2 п-1,2,3, 
+0.50[- W' (2IMAX-2IMIN-2+2JMAX-2JMIN-2) 
BE s qun м : 
Which was simplified to 
RESIDI = -(TERM)(W')+ U(non-island y') , (45) 
where 
TERM = IMAX - IMIN + JMAX - JMIN (46) 


The "hole relaxation" stipulated 2 individually 
there was no need for restriction on the values of island 
residuals, provided only that their sum, namely RESIDI, was 
eventually brought to zero. Therefore, when the absolute 
value of (45) was greater than the tolerable error € , a new 
value of the island ф', denoted pr", was defined by 

ое | (47) 
where CI was determined by requiring that the new island 
residual, RESIDI, be a negative fraction of RESIDI, that is, 


= 
RESIDI Ez -(aI) RESIDI , (48) 
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where QI was the appropriate over-relaxation coefficient for 
the chosen island. Since from (45) 

ши 0 (2 211 (05 
by substituting (47) into (49) and using (48), the constant 


CI was obtained to be 


(1.0 + al) 


sl TERM 


RESIDI (50) 


Ihe optimum value of the over relaxation coefficient, 
GI, was found to depend upon the size of the island selected 
for study. When the island sizes were small with respect to 
the scale of' the dominant wind stress, a value of 0.40 was 
assigned to al. If the island sizes were comparable to the 
scale of the stress, a smaller value of 0.20 was assigned. 
These time-tested selections were based upon the comparative 
need for more or less over-relaxation respectively as the 
island size was varied from small to large. 

More precisely, greater over-relaxation was neces- 
sary when the scale of the island was small in comparison 
to the scale of the wind stress and less over-relaxation was 
found necessary when the scale of the island was comparable 
to the scale of the wind stress. These results were to be 
expected since the scale of the motion, and hence the scale 
of the residuals will be determined by the scale of the wind 
stress forcing. 

Further details as to the application of this tech- 
nique in conjunction with that of the open oceanic areas 


can be found in Appendix A. 
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ТОО RESULTS 


А. А RESTATEMENT OF THE PURPOSE OF THE MODEL 

The main purpose of this model was to fully test and 
evaluate the boundary and relaxation techniques utilized. 
The accuracy of these results with the filtered model will 
be evaluated in later comparative studies with a free sur- 
face model. If the comparison is favorable, eventual incor- 
poration will be made in a multi-level world-ocean model 
with islands. 

The reader is reminded that the requirements of the 
model's solution are already required in more complicated 
multi-level baroclinic models. 

In the programming of the model, no special efforts were 
made to minimize computer time other than by selecting a 
value of the over-relaxation coefficient which was believed 
to be near the optimum value. With this in mind, the model 
required 28 minutes of C.P.U. time on the IBM 360/65 to 
integrate 364 time steps. Thus, with a time step of 
ДЭ? nate ly 14 hours, each of the 210 day integrations 
shown in the following section used less than one-half of 


one hour (.47 hours) of computer time. 


B. GENERAL DESCRIPTION OF ISLAND DATA CASE STUDIES 
Altogether, there were eight studies considered. Of 
these studies, there were seven different island positions 


examined, six of which were in the vicinity of the western 
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boundary; while the other island position was in the eastern 
boundary vicinity. The remaining study showed the island- 
less case. Because of the linearized form of the equations 
and the barotropic assumptions, the primary motivation for 
these eight studies was to specifically test and evaluate the 
model techniques and arrive at some quantitative idea as to 
what to expect in a more general, non-linear, multi-level 
future application of this model without excessive computer 
cost. 

Each island case study contained a variation of three or 
more island sizes centered about a specifically chosen grid 
location. For example, in the second data study, a specific 
grid position of "WESTERN GRID HALF-CENTERED" was chosen. 
The three distinct island sizes comprising this particular 
set were 400 kilometers square, 800 kilometers square, and 
1200 kilometers square. All three islands in this set were 
dat eredi in the western portion of the grid so that symmet- 
rical streamfunction patterns would appear to the north and 
south of the island. Each of the other island case studies 


was Similarly structured and can be referred to in Table l. 


С. EQUILIBRIUM CONSIDERATIONS AND FINDINGS 

A most important aspect of the numerical solution was to 
determine if the model had achieved a steady-state condition 
in each case study. Тһе ёо of an equilibrium condi- 
tion was accurately evaluated by means of a total kinetic 
energy СОА Са Е оп (Appendix A). At the end of each 


time step, this mathematical indicator projected an integral 
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picture of the steady-state characteristics of the evaluated 
island case. Specific illustrative comparisons between the 
steady-state characteristics of the islandless case and 
those of the island cases are shown in Figures 3 through 9. 
From an examination of these figures, it was clear that a 
statistical equilibrium was achieved in all cases by 210 days 
of integration. The streamline maps for the 210th day for 
each case is shown in Figures 10 through 34, while some of 
the data on each case is summarized in Table I. 

As shown in Table I, the following trends were noted: 
a) When the ‘island was located in a region of strong currents 
Шо о 574 6.7.8), the larger the size of the island, the 
lower the final value of the total kinetic energy.(Case 8-d 
was an unexplained exception to this rule.) When the island 
was located in a region of weak currents (Case 5), the strean- 
lines were easily bent around the island (Figure 23) pro- 
ducing a stronger flow in the central open ocean. Because 
the kinetic energy depends on the square of the current 
speed, this streamline bending by the island increased the 
total kinetic energy. b) Except for Case 8-d, the closer 
the island was positioned to the mainland in the westerly 
located studies, the smaller the final value of the total 
kinetic energy. dc) In Case 6 and in Case 7, the placement 
of a larger island in the path of a strong current reduced 
the final total kinetic energy value. 4) The comparative 
examination of Case 3 and Case 4 pointed out that the same 


final T.K.E. values were achieved by the symmetric placement 
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of equivalently sized islands south and north of the central 
latitude. In Case 4-b, an asymmetrically placed island 
Showed a different T.K.E. value which still conformed, how- 
ever, to those conditions stipulated above in (a) and (c). 
e) Finally, as might be anticipated, all island cases of 
comparatively small size had final T.K.E. values similar to 
that achieved in Case 1, the islandless case. 

As an additional guide in the above interpretation, the 
average kinetic energy was calculated from the T.K.E. by 
dividing by the number of non-island values which entered 
into the total. The results were unchanged; namely, that the 
larger the island, the smaller the average kinetic energy in 
all cases except the Eastern Island Case 5 and Case 8-b, 
where a large average kinetic energy appeared with the larger 
island. 

The contents of Table I display among other things the 


value of the streamfunction on the island (y. ) and maxi- 


island 
mum streamfunction value in the ocean basin Ор These 
values were significant in that they presented two quantita- 
tive measures of the current's volume transport per unit 
depth. Because of the zero normal flow boundary condition 


(14), ф 


) displayed the actual value of the volume trans- 
island 
port per unit depth between the island and the mainland. 
The reader is especially cautioned to note (in Case 8) that 


because of recirculation between the westerly boundary and 


the island, this value of а should be interpreted only 
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ав а measure of the net transport between an island and any 


mainland. 


The ae value accurately pointed out the maximum volume 


Cransport per unit depth to be found in the ocean basin. This 
value always was found in the area immediately adjacent to 
the western boundary of the mainland. 

As shown in Table I, in each data case study: a) The 


values of y. 


were smaller in i n 
713 г іп magnitude than those of 


В . b) When the size of the island was relatively larger, 


the magnitudes of y. 


island and prp were both smaller (Cases 


2, 3 and 4). c) When the island was positioned closer to 
the southern mainland boundary in the westerly located 


studies, the values of vy, 


9 and those of gre were both 


smaller. In the easterly located studies (Case 5), the 


Y 


values were nearly all the same, but the ША Ме values, 


island 
found near the western boundary of the ocean basin, consis- 
Benmtiy increased with the size of the island. This T 
increase seemed to be related to the increased gradient of 
the streamfunction near the central latitudinal axis of the 
basin; however, a further, more satisfactory explanation has 
not yet been found. d) In both Cases 3 and 4, the corres- 
ponding values of Vi sland and ‚жо were identical, again 
demonstrating the symmetrical characteristics of the linear 
solution. e) Finally, in Cases 6 and 7, smailer transports 
were recorded around the islands which were fully extended 
in the paths of strong currents, while the pe values 


found near the western ocean basin boundary were seen to be 
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Significantly smaller only in the case where the island was 
ШЫГУ extended north to south. 

Many features of the different solutions can only be 
seen in the streamfunction isoline maps which are displayed 
in Figures 10 through 34. 

Figure 10 showed the islandless case in which the numeri- 
cal solution compared favorably with the analytic solution of 
Munk. Figures 11 to 13 showed that when an island was placed 
in the center of an ocean gyre, the solution was not changed 
appreciably until the island size was about the same size or 
larger than the gyre itself. 

One of the interesting features found in Cases 3, 4 and 
ОШО ЕТ р итез 14 through 23) was that in this linear model, when 
an island was placed in the southern (northern) part of the 
basin, the latitude of maximum meridional transport shifted 
northward (southward). This shift was imperceptible in Case 
8 (Figures 30 through 34) in which the island had only a 
small north-south extent. There was some form of western 
intensification evident on the eastern shore of the island 
in every case study. Yet the maximum streamfunction values 
found on the eastern side of the island were always smaller 
in magnitude than the ае values depicted in Table I, 
which were always found adjacent to the western boundary of 
the basin. The island's eastern shore intensification was 
most evident in Case 6 (Figures 24 through 26) in which 


there was an anticyclonic gyre in the center of the basin, 


36 





Босго western intensified cyclonic gyres near the northern 


and southern ends of the fully extended island. 
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Y. CONCLUSIONS 


This numerical model, depicting the wind-driven circula- 
tion in an ocean basin containing an island, was successful 
in portraying the boundary conditions and relaxation tech- 
niques to be evaluated. Furthermore, a firm numerical еп 
has been established for a future non-linear 8 tee one and 
still later comparative study with a free surface model. The 
study of the combined relaxation technique employed in this 
model has established a much needed data resource which dis- 
plays the usefulness of the “Hole Relaxation" technique 


(Allen 1954) for the problem of the flow around an island. 
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